home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / tests / rtestode.mac < prev    next >
Encoding:
Text File  |  2003-02-09  |  6.3 KB  |  256 lines
  1. /* ODE tests */
  2.  
  3. /* Some of the tests are commented out.  The majority of them fail
  4.    due to some sort of testsuite interaction.  They pass if these
  5.    tests are run first, but then other tests fail */
  6.  
  7. kill(all);
  8. done;
  9.  
  10. /* Examples from "The Maxima Book" */
  11.  
  12. ode2(x^2*'diff(y,x)+3*x*y=sin(x)/x, y, x);
  13. y = (%C-cos(x))/x^3;
  14. ic1(%, x=1, y=1);
  15. y = -(cos(x)-cos(1)-1)/x^3;
  16. method;
  17. linear;
  18.  
  19. soln:ode2('diff(y,x,2) + y = 4*x, y, x);
  20. y = %K1*sin(x) + %K2*cos(x) + 4*x;
  21. method;
  22. variationofparameters;
  23. ic2(soln, x=0, y=1, diff(y,x)=3);
  24. y = -sin(x)+cos(x)+4*x;
  25. bc2(soln, x=0, y=3, x=2, y=1);
  26. y = -(3*cos(2)+7)*sin(x)/sin(2) + 3*cos(x) + 4*x;
  27.  
  28. ode2((3*x^2+4*x+2)=(2*y-1)*'diff(y,x), y, x);
  29. y^2-y = x^3+2*x^2+2*x+%C;
  30. method;
  31. separable;
  32.  
  33. /* Testsuite interaction causes this to fail */
  34. /* ode2(x^2*cos(x*y)*'diff(y,x) + (sin(x*y)+x*y*(cos(x*y)))=0, y, x);
  35. x*sin(x*y)=%C;
  36. method;
  37. exact; */
  38.  
  39. ode2( (2*x*y-exp(-2*y))*'diff(y,x)+y=0, y, x);
  40. x*exp(2*y) - log(y) = %C;
  41. method;
  42. exact; 
  43. intfactor;
  44. exp(2*y)/y;
  45.  
  46. /* testsuite interaction causes this to fail */
  47. /* ode2( 'diff(y,x)=(y/x)^2+2*(y/x), y, x);
  48. -(x*y+x^2)/y = %C;
  49. method;
  50. exact; */
  51.  
  52. ode2( 'diff(y,x)+(2/x)*y=(1/x^2)*y^3, y, x);
  53. y = 1/(sqrt( 2/(5*x^5) + %C)*x^2);
  54. method;
  55. bernoulli;
  56. odeindex;
  57. 3;
  58.  
  59. ode2( 'diff(y,x,2)-3*'diff(y,x)+2*y=0, y, x);
  60. y = %K1*exp(2*x) + %K2*exp(x);
  61. method;
  62. constcoeff;
  63.  
  64. ode2( 'diff(y,x,2)-4*'diff(y,x)+4*y=0, y, x);
  65. y = (%K2*x + %K1)*exp(2*x);
  66. method;
  67. constcoeff;
  68.  
  69. ode2(x^2*'diff(y,x,2)+x*'diff(y,x)-y=0, y, x);
  70. y=%K2*x-%K1/(2*x);
  71. method;
  72. exact;
  73.  
  74. ode2( x^2*'diff(y,x,2)+4*x*'diff(y,x)+2*y=0, y, x);
  75. y=%K1/x+%K2/x^2;
  76. method;
  77. exact; /*euler*/
  78.  
  79. ode2( x^2*'diff(y,x,2)+5*x*'diff(y,x)+4*y=0, y, x);
  80. y=(%K2*log(x)+%K1)/x^2;
  81. method;
  82. euler;
  83.  
  84. ode2( x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-1/4)*y=0, y, x);
  85. y=(%K1*sin(x)+%K2*cos(x))/sqrt(x);
  86. method;
  87. bessel;
  88.  
  89. ode2( x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-4)*y=0, y, x);
  90. y=%K1*%J[2](x)+%K2*%Y[2](x);
  91. method;
  92. bessel;
  93.  
  94. ode2( (x-1)^2*'diff(y,x,2)+(x-1)*'diff(y,x)+((x-1)^2-4)*y=0, y, x);
  95. y=%K1*%J[2](x-1)+%K2*%Y[2](x-1);
  96. method;
  97. bessel;
  98.  
  99. /* ode2( 'diff(y,x,2)+2*'diff(y,x)+y=exp(x), y, x);
  100. y=exp(x)/4+(%K2*x+%K1)*exp(-x);
  101. method;
  102. variationofparameters;
  103. /* Particular solution */
  104. yp;
  105. exp(x)/4; */
  106.  
  107. /* ode2( x*'diff(y,x,2)+('diff(y,x))^2=0, y, x);
  108. y='integrate(1/(log(x)+%K1),x)+%K2;
  109. method;
  110. freeofy; */
  111.  
  112. /* ode2( y*'diff(y,x,2)+('diff(y,x))^2=0, y, x);
  113. y^2/(2*%K1)=x+%K2;
  114. method;
  115. freeofx; */
  116.  
  117. eq: 'diff(y,x,2)+x*'diff(y,x)+exp(-x^2)*y=0;
  118. 'diff(y,x,2)+x*'diff(y,x,1)+%e^-x^2*y = 0;
  119. ans:ode2(eq,y,x);
  120. y = %k1*sin(sqrt(2)*sqrt(%pi)*erf(x/sqrt(2))/2)+%k2*cos(sqrt(2)*sqrt(%pi)*erf(x/sqrt(2))/2);
  121. is(ratsimp(ev(eq,ans,diff)));
  122. true;
  123. method;
  124. xformtoconstcoeff;
  125.  
  126. eq:x*'diff(y,x,2)+(x^2-1)*'diff(y,x,1)+x^3*y=0;
  127. x*'diff(y,x,2)+(x^2-1)*'diff(y,x,1)+x^3*y=0;
  128. ans:ode2(eq,y,x);
  129. y=%e^-(x^2/4)*(%k1*sin(sqrt(3)*x^2/4)+%k2*cos(sqrt(3)*x^2/4));
  130. is(ratsimp(ev(eq,ans,diff)));
  131. true;
  132. method;
  133. xformtoconstcoeff;
  134.  
  135. /* Tests of desolve */
  136.  
  137. eqn1:'diff(f(x),x) = sin(x)+'diff(g(x),x);
  138. 'diff(f(x),x,1) = 'diff(g(x),x,1)+sin(x);
  139. eqn2:'diff(g(x),x,2) = 'diff(f(x),x)-cos(x);
  140. 'diff(g(x),x,2) = 'diff(f(x),x,1)-cos(x);
  141. desolve([eqn1,eqn2],[f(x),g(x)]);
  142. [f(x)=%e^x*(at('diff(g(x),x,1),x = 0))-at('diff(g(x),x,1),x = 0)+f(0),g(x)=%e^x*(at('diff(g(x),x,1),x=0))-at('diff(g(x),x,1),x = 0)+cos(x)+g(0)-1];
  143. atvalue('diff(g(x),x),x = 0,a);
  144. a;
  145. atvalue(f(x),x = 0,1);
  146. 1;
  147. desolve([eqn1,eqn2],[f(x),g(x)]);
  148. [f(x) = a*%e^x-a+1,g(x) = cos(x)+a*%e^x-a+g(0)-1];
  149. remove(f,atvalue,g,atvalue);
  150. done;
  151.  
  152. atvalue('diff(g(x),x),x = 0,a);
  153. a;
  154. atvalue(f(x),x = 0,1);
  155. 1;
  156. desolve([eqn1,eqn2],[f(x),g(x)]);
  157. [f(x) = a*%e^x-a+1,g(x) = cos(x)+a*%e^x-a+g(0)-1];
  158.  
  159. /* eqn3: 'diff(f(x),x,2)+f(x)=2*x;
  160. 'diff(f(x),x,2)+f(x)=2*x;
  161. desolve(eqn3,f(x));
  162. f(x) = sin(x)*(at('diff(f(x),x,1),x = 0)-2)+f(0)*cos(x)+2*x; */
  163.  
  164.  
  165. /* Tests from:
  166.     F Postel and P Zimmermann, A Review of the ODE solvers of 
  167.     AXIOM, DERIVE, MACSYMA, MATHEMATICA, MUPAD and REDUCE
  168.     Procedings of the 5th Rhine Workshop on Computer Algebra
  169.     April 1-3, 1996, Saint-Louis, France
  170.     http://www.loria.fr/~zimmerma/ComputerAlgebra/
  171. */
  172.  
  173. /* Equation 1 - Linear polynomial */
  174. eqn: (x^4-x^3)*diff(u(x),x) + 2*x^4*u(x) = x^3/3 + C;
  175. (x^4-x^3)*'DIFF(u(x),x,1)+2*x^4*u(x) = x^3/3+C;
  176. ans:ode2(eqn,u(x),x);
  177. u(x) = ((2*x^3-3*x^2+6*C)*%E^(2*x)/(12*x^2)+%C)*%E^-(2*(LOG(x-1)+x));
  178. factor(ans);
  179. u(x)=%E^-(2*x)*(2*x^3*%E^(2*x)-3*x^2*%E^(2*x)+6*C*%E^(2*x)+12*%C*x^2)/(12*(x-1)^2*x^2);
  180. is(ratsimp(ev(eqn,ans,diff)));
  181. TRUE;
  182.  
  183. /* Equation (2) {firstord2} */
  184. eqn: -1/2*'diff(y,x) + y = sin(x);
  185. y-'diff(y,x,1)/2 = sin(x);
  186. ans: ode2(eqn,y,x);
  187. y = %e^(2*x)*(%c-2*%e^-(2*x)*(-2*sin(x)-cos(x))/5);
  188. ans: expand(ans);
  189. y = 4*sin(x)/5+2*cos(x)/5+%c*%e^(2*x);
  190. is(ratsimp(ev(eqn,ans,diff)));
  191. true;
  192.  
  193. /* Equation 3 - Linear Change of variables */
  194. eqn: 'diff(y,x,2)*(a*x+b)^2+4*'diff(y,x)*(a*x+b)*a+2*y*a^2=0;
  195. (a*x+b)^2*'diff(y,x,2)+4*a*(a*x+b)*'diff(y,x,1)+2*a^2*y = 0;
  196. ans:ode2(eqn,y,x);
  197. y = %k1*x/(a*x+b)^2+%k2/(a*x+b)^2;
  198. is(ratsimp(ev(eqn,%,diff)));
  199. true;
  200.  
  201. /* Equation 4 - Linear polynomial (Adjoint equation) 
  202.    See Zwillinger, p151
  203.    This crashed maxima-5.5
  204.    */
  205. eqn: (x^2-x)*'diff(y,x,2)+(2*x^2+4*x-3)*'diff(y,x)+8*x*y=1;
  206. (x^2-x)*'diff(y,x,2)+(2*x^2+4*x-3)*'diff(y,x,1)+8*x*y = 1;
  207. ans:ode2(eqn,y,x);
  208. false;
  209.  
  210. /* Equation 5 - Linear polynomial */
  211. eqn: (x^2-x)*'diff(y,x,2)+(1-2*x^2)*'diff(y,x)+(4*x-2)*y=0;
  212. (x^2-x)*'DIFF(y,x,2)+(1-2*x^2)*'DIFF(y,x,1)+(4*x-2)*y = 0;
  213. ans:ode2(eqn,y,x);
  214. false;
  215.  
  216. /* Equation 6 - Dependent variable missing */
  217. /* Another testsuite interaction problem */
  218. /* eqn:'diff(y,x,2)+2*x*'diff(y,x)=2*x;
  219. 'diff(y,x,2)+2*x*'diff(y,x,1) = 2*x;
  220. ans: ode2(eqn,y,x);
  221. y = sqrt(%pi)*%k1*erf(x)/2+x+%k2;
  222. is(ratsimp(ev(eqn,ans,diff)));
  223. true; */
  224.  
  225. /* Equation 7 - Liouvillian solution */
  226. eqn: (x^3/2-x^2)*'diff(y,x,2)+(s2*x^2-3*x+1)*'diff(y,x)+(x-1)*y=0;
  227. (x^3/2-x^2)*'diff(y,x,2)+(s2*x^2-3*x+1)*'diff(y,x,1)+(x-1)*y = 0;
  228. ode2(eqn,y,x);
  229. false;
  230.  
  231. /* Equation 8 - Reduction of order */
  232. eqn: 'diff(y,x,2)-2*x*'diff(y,x)+2*y=3;
  233. 'diff(y,x,2)-2*x*'diff(y,x)+2*y = 3;
  234. ode2(eqn,y,x);
  235. false;
  236.  
  237. /* Equation 9 - Integrating factors */
  238. eqn: sqrt(x)*'diff(y,x,2)+2*x*'diff(y,x)+3*y=0;
  239. sqrt(x)*'diff(y,x,2)+2*x*'diff(y,x)+3*y=0;
  240. ode2(eqn,y,x);
  241. false;
  242.  
  243. /* Equation 18 - Bernoulli equation */
  244. eqn: 'diff(y,x) + y = y^3*sin(x);
  245. 'DIFF(y,x,1)+y = SIN(x)*y^3;
  246. ans: ode2(eqn,y,x);
  247. y = %E^-x/SQRT(%C-2*%E^-(2*x)*(-2*SIN(x)-COS(x))/5);
  248. is(ratsimp(ev(eqn,ans,diff)));
  249. true;
  250.  
  251. /* Equation (20) {Clairaut} */
  252. eqn: (x^2-1)*'diff(y,x)^2 - 2*x*y*'diff(y,x) + y^2 -1 = 0;
  253. (x^2-1)*'diff(y,x)^2 - 2*x*y*'diff(y,x) + y^2 -1 = 0;
  254. ode2(eqn,y,x);
  255. false;
  256.